Growing Networks with Enhanced Resilience to Perturbation 



o 
o 

(N 

S 

in 



C 
i 

CO 



E 

c 

o 
o 



> 
l> 

O 
in 
o 

o 

-I— > 
o3 



c 

o 
o 



X 













+i 


+ 1 


-i 






Markus Bredcl and John Finnigan 
CSIRO Centre for Complex Systems, Canberra, ACT 2601, Australia 

(Dated: 5.5.2004) 

Scale- free (SF) networks and small world networks have been found to occur in very diverse contexts. 
It is this striking universality which makes one look for widely applicable mechanisms which lead 
to the formation of such networks. In this letter we propose a new mechanism for the construction 
of SF networks: Evolving networks as interaction networks of systems which are distinguished by 
their stability if perturbed out of equilibrium. Stability is measured by the largest real part of 
any eigenvalue of a matrix associated with the graph. We extend the model to weighted directed 
networks and report power law behaviour of the link strength distribution of the weighted graphs in 
the SF regime. The model we propose for the first time relates SF networks to stability properties 
of the underlying dynamical system. 



Recent studies have shown that a SF topology of the 
interaction network is a universal feature shared by many 
complex coupled systems. Examples are found in diverse 
fields including the WWW, traffic flow systems, social 
networks and genetic, metabolic, and protein folding net- 
works Q. Several mechanisms for the formation of SF 
networks are known. First SF networks can be built by 
preferential attachment 2] where new nodes form links 
preferentially to old nodes of high degree. Modifications 
to this procedure — e. g., incorporating an a priori as- 
signed individual node fitness — are still based on the 
general mechanism of peferential attachment and result 
in a slightly modified network topology |3|]. Further, SF 
networks can be considered as the result of a process op- 
timizing the diameter and link number in a network of 
given size Q , the direct construction of Hamiltonians for 
the network 5] or a thresholding mechanism 0. Very 
recently also, another mechanism via node merging was 
discovered pj. 

It has been shown that SF networks are very robust 
to random node removal thus allowing speculation 
that during evolution SF topologies might have been se- 
lected because of the inherent stability associated with 
their architectures. However, stability measures in the 
above work have been purely topological. The real situa- 
tion appears far more complicated. In most systems, the 
network topology only reflects the population dynamics 
as given by an underlying set of simultaneous equations. 
Heterogeneity in link strengths (denoting the strength of 
couplings in the equations) further complicates the pic- 
ture. 

Here, by directly relating a system's stability to its 
topology, we propose another mechanism to obtain SF 
networks. The mechanism we use comprises two essential 
steps: (i) random addition of new nodes to the network 
and (ii) selection of more stable networks, where stability 
is measured by the size of the largest eigenvalue of a 
matrix associated with the graphs. 

Consider a complex system the dynamics of which are 
described by some set of non-linear first order differential 
equations. If the system approaches an equilibrium state, 



FIG. 1: The attachment process. Dashed lines indicate nega- 
tive links, solid lines positive links. As described in the text a 
negative link has strength -1, a positive link strength +1. On 
the right hand side it is visualized how the new node links into 
the old network. It always forms one positive and one neg- 
ative in- and outlink with 4 randomly chosen vertices of the 
old network. The left hand figure shows how the attachment 
to the network changes the matrix M. 



a linear stability analysis can be undertaken, yielding 



Mx, 



(1) 



where x denotes deviations from equilibrium and M 
stands for the Jacobian matrix at equilibrium. The 
equilibrium is stable, if the largest real part A max of 
any eigenvalue of M is less than zero, where A max = 
max_\ go .( A f-)Re(A) and o~(M) is the spectrum of M. For 
generality we do not specify the exact form of these equa- 
tions underlying the dynamics of the network, but con- 
centrate on the Jacobian alone, following an approach 
similar to fgl. fToj| . 

The main diagonal elements of M arc set to ran — — F 
thus the populations are self regulated and normalized 
with respect to their intrinsic growth rates. Non-diagonal 
elements of M give the adjacency matrix of the network. 
We start allowing only matrix elements my £ { — 1, +1}, 
rriij = — 1 representing suppresion of node i's population 
by j and = +1 a stimulation to growth. 

To proceed, we propose a mechanism for the growth 
of a directed network with positive and negative connec- 
tions based on the A max < criterion. We start with 
a disconnected set of No — 4 nodes, hence m« = — 5a, 
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i,j = 1,...,4 and A max (4) = — 1, and continue with the 
iteration of the following steps: (i) Add a new node to the 
network, which forms I = 2 positive connections to ran- 
domly selected nodes, of which one is an in-link and the 
other an out-link. Next we add another I = 2 negative 
links in the same fashion [see Fig. (ii) The eigenval- 
ues of M are determined. Let A max (./V) denote the largest 
real part of all eigenvalues of M, and let A max (-/V — 1) be 
the same for the (N — 1) X (N — 1) matrix before insertion 
of the new node. If A max (A) > the configuration is re- 
jected immediately and we proceed with (i) and the last 
accepted configuration, (iii) If X max (N) < A max (A r — 1) 
the last node addition will be accepted. Otherwise, the 
acceptance probability is given by 

Accept = exp (-/3(A max (A) - A max (iV - 1))) , (2) 

where (3 is an inverse temperature-like parameter. Un- 
less the desired network size N has been reached, the 
algorithm continues with step (i). If before reaching the 
target network size a configuration to which no further 
node can be added is encountered, we start again from 
step (i). 

By construction every matrix M has the same number 
of positive and negative links. Furthermore, Tr(M) = 
£\Aj = —N giving A = 1/ATrAf = -1. Hence, in a 
stable matrix always Re(A) £ (0, — N) and, in particu- 
lar, Amax £ (0, —1]. Consequently, any broadening of the 
eigenvalue distribution in the direction of smaller eigen- 
values will also entail a larger eigenvalue closer to zero 
and thus a less stable system in our sense. Note, that by 
defining a Markov process in A max (iV) the above proce- 
dure also describes a Markov process in the space of all 
graphs with positive and negative links. The system size 
N gives the number of steps the walk has to perform. 

The algorithm can be interpreted in two ways. First, 
one may consider it only as a method to construct an en- 
semble of matrices (or graphs) with optimized stability 
properties. In this interpretation, the largest eigenvalue 
of a matrix might be considered as its 'energy', (3 being a 
measure for the fluctuations in the ensemble. Neglecting 
the system's growth, the procedure is then very similar 
to a Metropolis algorithm as commonly used for the nu- 
merical study of spin systems. 

On the other hand, it could also be conceived as a net- 
work evolution, mimicking a system's growth one node 
at a time. Every given interval of time a new 'species' is 
added to the system. This causes a perturbation to the 
population dynamics, that again settles into an equilib- 
rium. After relaxation, thermal fluctuations lead to per- 
turbations around the equilibrium which cause the less 
stable systems to collapse. In this view, the algorithm 
describes an evolutionary search, letting only the most 
stable systems survive. 

Degree distributions- In order to compare the results 
obtained by the algorithm to a random network evolu- 
tion, we calculate the degree distribution of the network 




in-degree 



FIG. 2: Example for the in-degree distribution of graphs (N = 
100) constructed with the above algorithm. The solid line 
(note the log scale on both axes) indicates a power law with 
exponent 7 « —2.23 ± 0.05. For comparison, the dashed line 
shows an exponential network as given by Eq. |QJ. The data 
are sampled for j3 = 100, iVo = 4, and represent averages over 
1000 independent runs. 



which would be obtained by only iterating step (i) of the 
algorithm. For simplicity we don't distiguish between 
positive and negative links in the evaluation of degree 
distributions in this letter. Following the rate equation 
approach one quickly obtains for nodes of degree 
larger than one (nodes of degree smaller than two remain 
from the initial conditions and don't affect the asymp- 
totic limit) 

{m(t + 1)> = (n<(t)) - (2/N(t)ni(t)) + (2/JV(i)n<_i(t)), 

(3) 

where Ui(t) stands for the number of nodes of degree 
i after iteration t. Assuming (rii(t)) ~ piN(t) in the 
limit of large network sizes, Eq.'s © yield an asymptotic 
degree distribution 

which is exponential. We now continue by comparing 
Eq. with degree distributions obtained in the ensem- 
ble of graphs constructed by the above algorithm. Figure 
12 shows simulation results for the in-degree distribution 
of networks of size 100 constructed with (3 — 100. Up 
to a finite size cut-off both the in- and out-degree distri- 
butions follow the same power law (rii) ~ i -7 with an 
exponent 7; n = 7 out = 7 = —2.23 ± 0.05. 
The influence of (3- Assuming a fixed network size the 
most important parameter in our construction is [3. To 
quantify the influence of (3 on the degree distribution we 
define a degree entropy Sdcg = J2iPi^°&Pi ( CI - -R- e ^ 0D- 

Figure illustrates changes in network structure with 
(3. For very small (J we find networks with an exponential 
degree distribution, which is well discribed by Eq. 
On the opposite end, for high (3, networks are SF. In the 
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FIG. 3: Simulation data showing changes in the degree en- 
tropy Sdeg with f3 for a system of size TV = 50. We find, that 
the relatively high degree entropies for small /3 correspond 
to exponential networks, while all networks constructed for 
high /3 are SF. In the intermediate region, a subensemble of 
the networks is SF, while some networks are neither SF nor 
exponential. 
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FIG. 4: Histogram counting the normalized frequency of how 
many times networks with A max occur if constructing 10 4 net- 
works of size 50. The distribution is bimodal with peaks at 
A = — 1 and A » —.65. We find that networks around the 
peak at A max = — 1 are SF. 



intermediate range, both SF and non-SF networks oc- 
cur. A more detailed investigation reveals, that the tran- 
sition between exponential and SF networks is related to 
a shift in typical largest real parts of eigenvalues. To 
illustrate this, Fig. 0] displays data for the distribution 
of A max . This distribution is bimodal. Corresponding 
to the two peaks we define two ensembles of networks, 
S- = {M e S|A max (M) < -A c }, A c = -.9, (the more 
stable subset) and S + = S — S- (the less stable sub- 
set). Interestingly, independent of 0, networks belonging 
to the more stable subset are SF while others belonging 
to S+ are typically not SF. A further finding is, that the 
degree distributions of the less stable networks are more 
likely to be exponential the farther away the second peak 
is from A = — 1. 

In the following we relate changes in j3 to changes in the 
weight of both peaks in the A max -distribution and thus 
to changes in network topology. Choosing a high value 
of f3, typical walks in A max are trapped in the vicinity 



of A = —1. Thus all networks are SF. As (3 is lowered, 
some walks escape the neighbourhood of A = —1 and 
give rise to a second peak in the A max -distribution. As 
(3 is further decreased more and more walks escape, thus 
shifting the balance towards the less stable non-SF net- 
works. Simultaneously, the smaller (3 the farther away 
from A = — 1 typical walks get. Finally, the choice of a 
very low (3 allows all walks to escape from A = — 1 and 
get close to A = 0, resulting in a dominating peak of less 
stable networks, which are exponential. 

Finding key mechanisms- It appears of interest to find 
which of the steps (i)-(iii) prove necessary to construct SF 
networks in the above way. For this we relax constraints 
in our construction procedure and study networks that 
result therefrom. 

We have already noted, that a high value of (3 is re- 
quired to obtain SF networks. Since, via Q, (3 deter- 
mines the chance that a positive step in the A-walks 
is accepted, a high (3 could also be interpreted as a 
trap confining walks to A = —1. Indeed, it turns out, 
that Eq. (J2J can be replaced by a sharp cut-off as 
Paccopt = 6 (A max (N) - Ai), A x < A c , still resulting in SF 
networks. One can conclude, that the functional form 
of Eq. (J2J is not essential. To construct 5_ it is only 
required to keep the eigenvalue distribution 'narrow' and 
'close' to A = — 1. 

As a next key element in our construction procedure 
we included random network growth. To check the im- 
portance of this step, we investigate the influence of op- 
timization for a narrow eigenvalue distribution alone. 
More specifically, we consider an uncorrelated random 
network of size N with L = 4(N — 1) links, of which 
half are negative and half positive. By swapping links 
between arbitrary nodes new network configurations are 
suggested, which are again rejected/accepted on the ba- 
sis of step (hi) and Eq. (J2J). Even choosing very small 
network sizes (N = 20) after 10 7 iterations we still find 
peaked degree distributions as in the initial random net- 
work. Hence, our simulations suggest that this procedure 
does not lead to the formation of SF networks. Thus, as 
in preferential attachment, growth seems to be a key re- 
quirement to obtain SF networks. 

Drawing link strength at random- So far we identified 
two necessary steps to produce SF graphs from the above 
matrix stability criterion: a tendency to trap eigenvalues 
close to A = — 1 and network growth. However, since we 
only allow links of strength +1 and —1 the procedure still 
lacks generality. In the following we draw link strengths 
from uniform distributions over the intervals [—1,0] and 
[0,1]. Hence we now have negative links of strength — 1 < 
s < and positive links in < s < 1, allowing for a 
discrepancy (mij)i^j ^ between the total weight of 
positive and negative links. In this way we construct an 
ensemble of weighted directed graphs 

As for the graphs with link strength s = ±1 we again 
find a bimodal distribution of A max (iV). Choosing a suit- 
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FIG. 5: Distribution of the strengths of links for the more 
stable subset S- for networks of size N — 50 constructed 
with /3 = 50 (points). Interestingly, the more stable graphs 
have more strong links than the less stable ones (not shown) . 
In the ensemble of more stable graphs the distribution follows 
two distinct power laws Pr(s) ~ s~ s ~/+ (dashed lines), with 
exponents 8+ = .51 ±0.03 (s > 0) and 5- = .45 ±0.02 (s < 0). 
The inlet shows the different power laws holding for s > and 
s < on log-log scales. The data are logarithmically binned 
and represent averages over 10 independent runs. 



able A c a more stable (S_) and a less stable ensemble 
(S+) of graphs can be defined as previously. Again, net- 
works belonging to the more stable ensemble are SF (with 
exponents y ln = j out = —2.35 ± 0.04 for N — 100 and 
(3 = 50). Similarly, graphs in the less stable ensemble 
approach an exponential degree distribution for small (3. 

To proceed, we examine the distribution of link 
strength of the more stable graphs S- [see Fig. 
Generally, as would be expected, weak links are much 
more frequent than strong links. Further, in the case of 
the ensemble of more stable matrices one finds distinct 
power laws Pr(s) ~ s~ s -/+ for both the negative and 
the positive branch of the distribution. For s < an 
exponent <5 + = .51 ± 0.03 is obtained, while for s < it 
holds 6- — .45 ± 0.02. A similar behaviour, i.e. domi- 
nance of relatively weak links and the existence of only a 
few strong links ('hot-lin ks') , is expected in many empir- 
ical networks in biology [l5j. Further, a recent empirical 
study confirms power law behaviour in the link strength 
distribution of some SF networks [l4|. 

Most of the empirically investigated networks are 
found to be not only SF, but to also show a high amount 
of 'cliquishness' pj. This is measured by c, the clustering 
coefficient, which calculates the fraction of links between 
neighbours' neighbours which are actually present. In 
order to factor out the effects of a specific degree distri- 
bution, the values (c) found in the ensemble S- are com- 
pared to averages over an ensemble of randomized graphs 
(crand) with the same degree distribution (cf. Il2j) and 
to random networks of the Erdos-Renyi type [l3j (cer)- 
For (3 = 50, N = 100 and link strength |s| randomly 
drawn from [0, 1] our simulations result in (c) = .078 
to be compared with (c ran d) = -045 and cer = .025. 



Hence, networks in the stable ensemble are substantially 
more cliquish than random networks. Preliminary stud- 
ies seem to indicate only a slight dependence on (3 and 
a ratio (c)/(c ranc i) which grows with the system size N. 
Similar experiments for the average shortest path length 
give values very close to that found in random networks, 
hinting to a small world like topology of our graphs. 

In summary, using a stability criterion based on the 
largest real part of eigenvalues of a matrix associated 
with a graph, we have presented a model to generate an 
ensemble of graphs distinguished by their stability. In 
this context, more stable graphs turn out to have SF 
degree distributions. We identified two key mechanisms 
to obtain SF graphs in this way: growth and a tenden- 
dency to keep the real part of the eigenvalue distribu- 
tion 'narrow'. Extending the model to generate weighted 
graphs, we find stable graphs that are SF and exhibit a 
power law distribution of the link strengths. We find that 
the networks are substantially more cliquish than random 
networks, while still exhibiting a path length very sim- 
ilar to random networks. The model thus for the first 
time relates SF 'small world' like graphs to a stability 
criterion which is directly associated with the underlying 
equations' dynamics and thus adds another notion to the 
understanding of SF networks as networks distinguished 
by robustness against perturbations. 
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